Atomistic simulations of structural and thermodynamic properties of bilayer graphene 
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We study the structural and thermodynamic properties of bilayer graphene, a prototype two- 
layer membrane, by means of Monte Carlo simulations based on the empirical bond order potential 
LCBOPII. We present the temperature dependence of lattice parameter, bending rigidity and high 
temperature heat capacity as well as the correlation function of out-of-plane atomic displacements. 
The thermal expansion coefficient changes sign from negative to positive above w 400 K, which is 
lower than previously found for single layer graphene and close to the experimental value of bulk 
graphite. The bending rigidity is twice as large than for single layer graphene, making the out-of- 
plane fluctuations smaller. The crossover from correlated to uncorrelated out-of-plane fluctuations 
of the two carbon planes occurs for wavevectors shorter than « 3 nm -1 . 

PACS numbers: 81.05.ue, 61.48.Gh, 65.80.Ck, 87.16.D- 



I. INTRODUCTION 



Bilayer {BL) graphene has unique electronic proper- 
ties and its chiral quasiparticles with parabolic dispersion 
make it different from both single layer (SL) graphene 
and bulk graphite 1 . The energy gap of BL graphene 
can be opened and tuned by applying a voltage, with 
promises for applications 2,3 . Also the possibility of some 
exotic many body phenomena, such as pseudospin mag- 
netism^, have been discussed. For these reasons, BL 
graphene is currently subject of great interest. However 
the knowledge of its structural properties is still very 
poor. It was shown experimentally that, BL graphene 
is also corrugated 5 like SL graphene, but no systematic 
study has been carried out. This corrugation (ripples) 
may constitute an important scattering mechanism for 
electrons^ and ripples can give rise to charge inhomo- 
geneities (electron and hole puddles )£. Although impor- 
tant for their relation to electronic properties, the struc- 
tural properties of BL graphene are also important from 
the point of view of statistical mechanics since the BL 
graphene is a unique realization of crystalline membranes 
formed by two atomic layers. 

Assessing the structure of a BL graphene is experi- 
mentally challenging and theoretical calculations can be 
particularly helpful. Since the observed corrugations are 
on a scale much larger than interatomic distances, ab- 
initio simulations are not feasible. This interesting range 
of lengths (e.g. for electron interactions with ripples) 
is, however, not necessarily well described by continuum 
medium theories 8 . Atomistic simulations based on accu- 
rate empirical interaction potentials are particularly suit- 
able for this purpose. We have recently studied the struc- 
tural and thermodynamic properties of SL graphene^— 
by Monte Carlo (MC) simulations based on the LCBOPII 
bond order potential^ 2 -. Here we present the results of 
similar calculations for BL graphene, where a new as- 
pect related to the correlation of atomic displacements 
in different layers arises. 



II. METHOD OF CALCULATION 



We perform MC simulations in the NPT ensemble 
at pressure P = and temperature T with periodic 
boundary conditions for samples of N = 16128 and 
N = 8640 atoms per layer. When not specified, the 
presented results are for the largest sample. The equi- 
librium size at T = K of the N = 16128 sample is 
L x = 20.66 nm in x and L y — 20.448 nm in y direction 
and that of the N = 8640 sample is L x = 14.757 nm and 
L y = 15.336 nm. The finite size of our sample defines 
the lowest accessible wavevectors is x and y directions as 
q x = 2n/L x and q y — 2n/L y . Motivated by the results 
of recent quantum MC calculations^, we have slightly 
modified the long-range part of LCBOPII as to have an 
interlayer binding energy of 50 meV/atom a gain st the 
25 meV/atom of the parametrization of Ref. |l2|, while 
keeping the interlayer compressibility constant. 

We equilibrate the sample for at least 5 ■ 10 5 steps 
(1 MC step corresponds to N attempts to a coordi- 
nate change), using the recently introduced MC sampling 
based on collective atomic moves (wave moves)ii in addi- 
tion to conventional MC moves. This technique was suc- 
cessfully introduced for SL graphene. For BL graphene 
it was extended as follows. Wave moves are applied to 
both layers simultaneously, or only to the upper or lower 
layer, with equal probabilities for the three cases. The 
amplitude A\ of the wave moves applied to both layers 
simultaneously is different from the amplitude A2 of the 
wave moves applied to either upper or lower layer sepa- 
rately. The amplitudes A\ and A2 are chosen in such a 
way that the acceptance rate for wave moves is between 
0.4 and 0.5 for any of these three cases. 

Further 5 • 10 5 MC steps are used to evaluate the tem- 
perature dependence of the ensemble averages. 
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FIG. 1. (color online) Temperature dependence of the in- 
plane lattice parameter a of SL (circles, solid, blue, from 
Ref. Hoi ) and BL (circles, dashed, red) graphene, and of the 
interlayer distance c of BL graphene (diamonds, dash-dotted, 
green). At T = 0, a SL = 0.24595 nm, a BL = 0.24583 nm, 
c = 0.33371 nm. 



III. RESULTS 

The temperature dependence of the in-plane lattice pa- 
rameter a and of the interlayer distance c of BL graphene 
are shown in Fig. [TJ The in-plane lattice parameter a of 
BL graphene decreases with increasing temperature up to 
about 400 K, yielding a negative thermal expansion co- 
efficient a a = d\na/dT = (-3.0 ± 0.7) • 10~ 6 K" 1 in the 
range 0-300 K. The behavior of a(T) differs from that of 
SL graphene, which has a minimum of a at T ps 900 K 
and a a = (-4.8 ± 1.0) • 10~ 6 KT 1 (see Ref. Q3) in the 
range 0-300 K, and is similar to bulk graphite, which 
has a minimum of a between 300 and 500 Ki£>±&. We 
note that our approach is classical and therefore not ap- 
propriate in the low temperature limit. However, since 
the thermal expansion is mostly determined by the low- 
frequency bending modest, a classical description is al- 
ready justified below room temperature. Indeed for single 
layer graphene, our results for a(T) between 100 K and 
400 K agree very well with those of Ref. [H] where the 
quantum statistics of phonons was taken into account. 

In Ref. [TJ, the temperature dependence of a for SL 
graphene and bulk graphite has been determined in the 
quasiharmonic approximation with phonon frequencies 
and Gruneisen parameters calculated from first princi- 
ples. While for the case of bulk graphite these calcu- 
lations reproduce the non monotonic behavior of a(T) 
observed experimentally, for SL graphene a(T) keeps de- 
creasing up to high temperatures. In our simulations of 
SL graphene 1 ^ we found instead a non monotonic be- 
havior of a(T). The experimental value of a(T) for SL 
graphene that was measured up to 400 K— seems to sup- 
port our results. 

The discrepancy with quasiharmonic results should 
be due to the fact that this method 1 ^ neglects self- 




FIG. 2. Schematic view of BL graphene (solid lines), hi 
and hi are out-of-plane deviations with respect to the middle 
planes (dashed lines). The unit vectors n\ and ni are the nor- 
mals to each point in the upper and lower layer respectively. 
no is the normal to the reference plane, c is the interlayer 
distance. The figure is schematic and does not show the real 
scale of the fluctuations. 

anharmonic effects 1 ^, namely multiphonon contributions 
to the free energy. Of course, in our simulations, the 
thermal expansion is calculated directly and all anhar- 
monic effects are taken into account. Unfortunately, we 
do not have results for bulk graphite with the same in- 
plane area, due to the long range part of our potential 
that, with a cut off of 0.6 nm, requires to simulate sam- 
ples with at least four layers with periodic boundary con- 
ditions. Nevertheless, we believe that the fact that the 
thermal expansion of BL graphene is similar to the one 
resulting from quasiharmonic theory for bulk graphite 
suggests that multiphonon processes are much less im- 
portant in BL graphene, compared to SL graphene. 

In Fig. [TJ we also show the interlayer distance c that 
grows with temperature, similarly to bulk graphite 14 , 
with an out-of-plane thermal expansion coefficient a c = 
d\nc/dT = (3.5 ± 0.5) • 10" 5 K"\ which is compara- 
ble to the experimental value for bulk graphite, a c = 
2.7- 10~ 5 K- 1 (see Ref.[H). 

We now proceed to a study of thermal bending fluc- 
tuations. In the continuum limit graphene can be de- 
scribed as a flexible crystalline membrane&^iii which is 
characterized by a two component in-plane phonon field 
U a (x),a — 1,2 and a one component out-of-plane dis- 
placement field h(x). The effective free energy is given 
by the sum of bending energy and in-plane elastic energy^ 




where the strain tensor u a p is 

u a p = ^ (daup + dpu a + dji dph) , (2) 

k is the bending rigidity and (i and A are Lame coeffi- 
cients. 
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In first approximation, BL graphene can be consid- 
ered as two SL graphene layers interacting with each 
other. The natural way to describe BL graphene, is to 
use the out-of-plane deviations from the center of mass 
of each layer, h± and /12 for upper and lower layer re- 
spectively as sketched in Fig. [2j Thus, BL graphene 
can be parametrised by the average height fluctuation 
field h = (hi + h<z) /2 and thickness fluctuation field 
Sh = hi — hi- 

The part of the Hamiltonian (p} related to out-of-plane 
displacements can thus be written as: 



Hout — 



l 2 x (^k ( 



V 2 V 



K (V 2 /l 2 



27 (Shy 



(3) 

where the first two terms are responsible for the bending 
energy of the upper and lower layers and K is the bend- 
ing rigidity per layer. We have introduced the last term 
characterized by the parameter 7 to account for inter- 
layer interactions. Substituting hi and /12 with /i±<5/i/2, 
we obtain: 
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In the harmonic approximation, which means neglect- 
ing the last term of the strain tensor $2$, the out-of-plane 
h(x) and in-plane u a (x) modes are decoupled. In this ap- 
proximation, the mean square Fourier components of the 
field h(q) with wavevectors q are: 



and of the field 5h(q) are: 
<l*%1| 2 > = 



TV T 
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(5) 



(6) 



where N is the number of atoms per layer and Sq is the 
area per atom in the layer. If the bending rigidity of a SL 
graphene is the same as the bending rigidity per layer of 
BL graphene, then it follows from Eq. (O, that (\h(q)\ 2 ) 
for BL graphene is twice smaller than for SL graphene. 
This is actually a very good approximation as we will 
show below. 

We further introduce the notation H(q) = (\h(q)\ 2 ) 
and AH(q) = (\Sh(q)\ 2 ). 

An alternative way to describe out-of-plane fluctua- 
tions is via the unit vector normal to the average surface 
between two layers: 



rii(x) 



dih 
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(7) 



with i = 1,22. 

The correlation function of the normals, G(q) = 
(\n(q)\ 2 ) is equal to q 2 H(q) if \Vh\ 2 < 1. Thus, in the 
harmonic approximation 
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FIG. 3. (color online) Normal-normal correlation function 
G(q)/N at T = 300 K for SL (solid blue line) and BL (dashed 
red line) graphene compared to q 2 H(q) for single (dotted ma- 
genta line) and BL (dash-dotted red line) graphene. The solid 
black straight line shows the fit (~ q~ 2 ) in the harmonic part. 



which is a factor 2 smaller than G(q) in SL graphen o 9 ' 11 . 

The correlation functions H(q) and G(q) are calculated 
independently as described below. In principle to calcu- 
late H(q), we have to calculate the Fourier transforms 
of the atomic displacements h(x). However, the atomic 
positions in a generic configuration in MC simulations 
are discontinuous and should be smoothed. This prob- 
lem is related to the numerical calculations of derivatives 
and different operators on the hexagonal lattice^. Our 
procedure is the following. Let ho be the z-coordinate 
of an atom and h a , h^ and h c the z-coordinates of its 
three nearest neighbors. Then the averaged out-of-plane 
displacement of the central atom ho is defined as: 



h 



ho + \ (h a + h + h c ) 



(9) 



This value is used to calculate the Fourier components 
h(q) using the wave vectors defined by periodic boundary 
conditions of the undistorted lattice. The normals needed 
to calculate G(q), instead, are automatically smooth be- 
cause they are calculated as averages of the normals to 
the three planes defined by three vectors, connecting the 
central atom to its three nearest neighbors 9 -^. For BL 
graphene, we calculate the correlation function G(q) for 
the normals of all atoms in the two layers. 

Fig. |3] shows the correlation functions G(q)/N and 
q 2 H(q)/N for SL and BL graphene at T = 300 K. We 
plot these functions as a function of q = \q\ by giv- 
ing their average value at all allowed wavevectors with 
the same modulus. The difference between G(q)/N and 
q 2 H(q)/N is negligible for q < 10 nrn -1 where the con- 
dition \Vh\ 2 <gc 1 is satisfied. The functions H(q) and 
G(q) behave according to the harmonic approximation 
Eqs. ([S]) and © for q from 3 to 9 nrn -1 as it is also shown 



4 




-l 

q, nm 



FIG. 4. (color online) Average height (|/t(g)| 2 ) (blue lines) and 
thickness {\Sh(q)\ 2 } (red lines) fluctuations of BL graphene at 
T = 300 K (dashed lines) and T = 1500 K (solid lines). Black 
solid lines show the fit according to the Eqs. (O-©. 

in Fig.[3l In this interval the correlation functions for BL 
graphene are about twice smaller than for SL graphene, 
which means that the effective bending rigidity for BL 
graphene is twice larger than the one of SL graphene, as 
we had guessed above. The deviation from the harmonic 
approximation for q < 3 nm -1 is due to the coupling 
between bending and stretching modes in Eq. @ 8 . 

Fig. @] shows the correlation functions H(q)/N and 
AH(q) /N of BL graphene for T — 300 K and T = 1500 K 
together with the harmonic fit according to Eqs. (H])-©. 
The function AH(q) is specific of a bilayer and has no 
analog for single layer membranes. First of all, we note 
that AH(q) follows the harmonic approximation Eq. © 
in the whole studied range of q, even where the devia- 
tions of H(q) from the harmonic approximation Eq. ([5]) 
are pronounced. This means that thickness fluctuations 
are much less coupled to in-plane fluctuations, than av- 
erage out-of-plane fluctuations. 

The second noticeable point, is the g-independent be- 
havior of AH(q) for q < q* rj 3 nm . In turn, this 
means that, in this range of q, the out-of-plane fluctua- 
tions of the two carbon layers are strongly coupled and 
only one soft mode h(q) survives. Therefore, at scales 
larger than 2ir/q* sa 2 nm, BL graphene can be con- 
sidered as a single membrane, whereas at smaller scales, 
hi(q) and h2(q) fluctuate rather independently. Indeed 
it follows from Eqs. ([5]), © that if one neglects the in- 
terlayer coupling 7 in Eq. (|6]) one has: 

N T 

(M?)M<z)> = 0, <M?)| 2 > = (\h 2 (q)\ 2 ) = — (10) 

Sq nq 4 

In general, the perfect coincidence of AH(q) calculated 
from the MC simulations, with the theoretical predic- 
tion ^ of the Hamiltonian (Q} confirms the correct choice 
of Hamiltonian to describe BL graphene. 
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FIG. 5. (color online) Temperature dependence of the bend- 
ing rigidity « of SL graphene (circles, solid, blue), bending 
rigidity per layer k of BL graphene (circles, dashed, red) and 
parameter 7 of BL graphene (diamonds, dash-dotted, green). 



The crossover at q* from independent to coherent fluc- 
tuations in the two layers is important for the scatter- 
ing of electrons from height fluctuations in BL graphene, 
which is determined mainly by long range fluctuations 
with strongly (independent correlation functions (com- 
pare with Ref. [6| for SL graphene). Therefore, fluctu- 
ations of the interlayer distance become irrelevant for 
electrons with wavevector k < q*. Moreover, for sam- 
ple sizes L > 2n/q* w 2 nm the height fluctuations in BL 
graphene are expected to be weaker than in SL graphene, 
because in the regime of coherent fluctuations the bilayer 
is twice stiffer than a single layer. This is qualitatively 
confirmed by the results presented in Fig. [5] where we 
compare the values of (h 2 ) for SL and BL graphene. 

The temperature dependence of the parameters k and 
7 of BL graphene are presented in Fig. [5] together with 
the parameter k of SL graphene. The parameter 7 of BL 
graphene decreases with temperature, which is not sur- 
prising. This parameter is responsible for the interlayer 
coupling, and it decreases with temperature since the in- 
terlayer distance c increases with temperature (Fig. [TJ . 
The effective bending rigidity k grows with temperature 
in agreement with the general theory of crystalline mem- 
branes^, as well as with our previous numerical results 
for SL graphen o 9 ! 11 . The behavior of liquid membranes 
is known to be opposite, with k decreasing with temper- 
ature^. The statement that k decreases with T also for 
graphene^i is therefore in disagreement with general ar- 
guments&i 9 - and our results. The point is that the origin 
of the main anharmonic effects in liquid and crystalline 
membranes are completely different. For liquid mem- 
branes they originate from high order terms of the mean 
curvature in V/i, which results in perturbative correc- 
tions to k that arc of the form Thiqa < with a the 
interatomic distanc e 19 ! 20 . For crystalline membranes, in- 
stead, perturbative corrections to k(T) due to the cou- 
pling of bending and out-of-plane fluctuations are much 
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FIG. 6. (color online) Height fluctuations of SL graphene 
(solid blue line) compared to one of BL graphene (dashed red 
line) as a function of MC step at T = 300 K. 
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FIG. 7. (color online) Temperature dependence of the molar 
heat capacity at constant volume Cv of SL (solid blue line) 
and BL (dashed red line) graphene. Data obtained for N = 
8640 atoms sample. 



stronger, positive and proportional to T/q 2 (Ref. 19). 

Actually the fact that dn/dT > for crystalline mem- 
branes has a very simple meaning: as the temperature 
increases, the amplitude of corrugation also increases, re- 
sulting in a strengthening of the membrane^. As already 
mentioned, the bending rigidity per layer of BL graphene 
turns out to be very close to that of SL graphene, which 
is not surprising since the interlayer coupling is much 
weaker than the in-plane chemical bonding. However, 
since the renormalization of n is strongly g-dependent 
for crystalline membranes, the definition of k(T) should 
be further specified. What is shown as k(T) in Fig. [5j 
and what was previously calculated for SL graphene in 
Refs^iii are the results of a best fit of the correlation 
functions G(q) and H(q) in the q range where the slope 
can be well approximated by the harmonic behavior of 
Eqs. (U]) and (|5J|. Since, in this interval of q, the out-of- 
plane fluctuations of either layer of BL graphene are of 
the same order as those of SL graphene (see Eq. (jTU)) ). it 
is not surprising that the temperature dependence of k 
for BL graphene is only marginally smaller than for the 
one of SL graphene. 

It is important to notice, however, that the macro- 
scopic behavior of the bending rigidity of free membranes 
for q — >• at finite temperature is divergent as q~ v with 
j] w 0.85 (see Refs. H, [ill). The size of the BL graphene 
samples used here makes an estimate of r\ for this case not 
precise enough as to be compared quantitatively to that 
found for SL graphene^i, but the qualitative behavior 
shown in Fig. [3] is very similar for SL and BL graphene. 

The mean square height fluctuations (h 2 ) — ^ q H(q) 
are equally size-dependent. Since the sum over q is diver- 
gent at the lower limit q m in — 27r/ L, (h 2 ) is mostly deter- 
mined by the effective n(q) for the smallest wavevectors 
and therefore, for large enough samples, it should scale 
as L 2 ~ v (see Refs.l^.lllr). According to Fig.|U deviations 
of H(q) from harmonic behavior occur for q < 1 nm _1 



and thus, the crossover from harmonic behavior h 2 oc L 2 
to the anharmonic one h 2 oc L 2 ~ v takes place for sample 
size L ks 6 nm. 

To characterize qualitatively the anharmonicity at the 
atomic scale, we calculate the temperature dependence 
of the molar heat capacity at constant volume 



3i? 

IT 



dU 



(ii) 



where U is the potential energy and R the gas con- 
stant. In Fig. [7] we compare the results with those of 
SL graphene^. 

Three and four phonon processes result in the lin- 
ear growth of Cy at high temperatures 15 . One can see 
that SL and BL graphene are almost the same as ex- 
pected, since phonons of the whole Brillouin zone con- 
tribute to this quantity and the phonon spectra of SL 
and BL graphene differ only slightly close to the T point 
(see, e.g., the calculated phonon spectra of graphene and 
graphite in Ref. [H]) . 



IV. SUMMARY 

In conclusions, we have studied several temperature 
dependent properties of BL graphene by means of classi- 
cal MC simulations. The high temperature heat capacity 
is similar to that of SL graphene, whereas the thermal 
expansion is essentially different and close to the one ex- 
perimentally observed in graphite. 

We also introduced a new Hamiltonian which accounts 
for interlayer interactions in BL graphene and showed 
that it correctly describes the behavior of BL graphene. 
We have found that, depending on the wavevector, the 
height fluctuations in the two layers are either coher- 
ent (for q < q*) or incoherent (for q > q*) with q* w 
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3 nm" 1 at room temperature and we have discussed the 
consequences of this fact for observable properties, like 
height fluctuations and electron scattering. 
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